Calculate predation and spatial overlap with prey

Load packages & source functions

library(sdmTMB)
library(tidyverse)
library(tidyverse)
library(viridis)
library(devtools)
library(terra)
library(tidylog)
library(RColorBrewer)
library(patchwork)
library(egg)
library(ggpmisc)
library(ggh4x)
library(ggtext)

# Set path
home <- here::here()

# Source functions for overlap, predation and plotting
for(fun in list.files(paste0(home, "/R/functions"))){
  
  source(paste(home, "R/functions", fun, sep = "/"))
  
}

# Palette for plotting
pal <- brewer.pal(name = "Paired", n = 8)

Load prediction grid

# scale environmental variables with respect to catch data, and diet covariates with respect to diet data
catch <- read_csv(paste0(home, "/data/clean/catch_clean.csv")) |> drop_na(depth, oxy, salinity, temp)
diet <- read_csv(paste0(home, "/data/clean/stomachs.csv"))

pred_grid <- 
  bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
            read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv"))) |> 
  filter(quarter == 4) |> # Not needed in theory for saduria...
  drop_na(saduria) |> 
  mutate(year_f = as.factor(year),
         quarter_f = as.factor(quarter),
         month_f = as.factor(month),
         oxy_sc = (oxy - mean(catch$oxy)) / sd(catch$oxy),
         temp_sc = (temp - mean(catch$temp)) / sd(catch$temp),
         temp_sq = temp_sc^2,
         depth_sc = (depth - mean(catch$depth)) / sd(catch$depth),
         depth_sq = depth_sc^2,
         salinity_sc = (salinity - mean(catch$salinity)) / sd(catch$salinity),
         pred_length_sc = 0,
         saduria_sc = (saduria - mean(diet$saduria))/sd(diet$saduria))

Load models

# Density
mcod <- readRDS(paste0(home, "/output/mcod.rds"))

# Diet
mher <- readRDS(paste0(home, "/output/mher.rds"))
mspr <- readRDS(paste0(home, "/output/mspr.rds"))
msad <- readRDS(paste0(home, "/output/msad.rds"))

Predict with sims from density and diet models to propagate uncertainty

# predict and then for loop to summarise and avoid having too large pred grids in the memory
nsim <- 500

sim_dens <- predict(mcod, newdata = pred_grid, nsim = nsim) |> as.data.frame()

sim_her <- predict(mher, newdata = pred_grid, nsim = nsim) |> as.data.frame()
sim_spr <- predict(mspr, newdata = pred_grid, nsim = nsim) |> as.data.frame()
sim_sad <- predict(msad, newdata = pred_grid, nsim = nsim) |> as.data.frame()

The density prediction is used later for spatial overlap + it gets left_joined to the other sims, so it gets a special treatment

pred_grid_density <- pred_grid |> 
  bind_cols(sim_dens) |> 
  pivot_longer(cols = starts_with("V"), names_to = "sim", values_to = "sim_density") |> 
  dplyr::select(X, Y, year, saduria, biomass_spr, biomass_her, sim_density, sim) |> 
  mutate(sim_density = exp(sim_density)) |> 
  filter(sim_density < quantile(sim_density, 0.9999))

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, sim_density) [was 504897x532, now 252448500x34]

mutate: changed 252,448,500 values (100%) of 'sim_density' (0 new NA)

filter: removed 25,245 rows (<1%), 252,423,255 rows remaining

Plot annual predation with uncertainty

# A prediction grid cell is 3*3 km
area <- 9

# Summarise predator biomass by year and sim and join that to calculate per capita predation
pred_dens <- pred_grid_density |> 
  summarise(cod_biomass = sum(sim_density*area),
            .by = c(year, sim)) 

# Apply functions and combine
tot_pred <- bind_rows(
  get_annual_predation(sim_her, threshold = 1, power = 2, prey_name = "Herring"),
  get_annual_predation(sim_sad, threshold = 1, power = 2, prey_name = "Saduria"),
  get_annual_predation(sim_spr, threshold = 1, power = 3, prey_name = "Sprat")
)

# Prepare data for fitting gams and for plotting
clean_pred <- tot_pred |>
  pivot_longer(c(-year, -species)) |>
  mutate(var = ifelse(grepl("_sd", name), "sd", "mean"),
         metric = ifelse(grepl("cap", name), "Per capita (kg/kg)", "Total (tonnes)")) |>
  dplyr::select(-name) |>
  pivot_wider(values_from = value, names_from = var) |>
  mutate(species = as.factor(species))

# Fit gamma-GAMs to annual predation metrics
gamma_gam_tot <- sdmTMB(mean ~ species + s(year, by = species),
                        spatial = "off", 
                        family = Gamma(link = "log"),
                        data = clean_pred |>
                          filter(metric == "Total (tonnes)") |> 
                          mutate(mean = mean/1000))

gamma_gam_cap <- sdmTMB(mean ~ species + s(year, by = species),
                        spatial = "off", 
                        family = Gamma(link = "log"),
                        data = clean_pred |>
                          filter(metric == "Per capita (kg/kg)"))

nd <- expand_grid(species = as.factor(c("Sprat", "Herring", "Saduria")),
                  year = unique(clean_pred$year))

p_tot <- predict(gamma_gam_tot, newdata = nd, se_fit = TRUE)
p_cap <- predict(gamma_gam_cap, newdata = nd, se_fit = TRUE)

p <- bind_rows(p_tot |> mutate(metric = "Total (tonnes)"),
               p_cap |> mutate(metric = "Per capita (kg/kg)"))

clean_pred <- clean_pred |> 
  mutate(lwr = (mean - sd),
         upr = (mean + sd),
         lwr = ifelse(lwr < 0, 0, lwr),
         mean = ifelse(metric == "Total (tonnes)", mean / 1000, mean),
         lwr = ifelse(metric == "Total (tonnes)", lwr / 1000, lwr),
         upr = ifelse(metric == "Total (tonnes)", upr / 1000, upr))

p0 <- ggplot(clean_pred, aes(year, mean)) + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0, color = "grey40", alpha = 0.7,
                linewidth = 0.4) + 
  geom_point(color = "grey30", alpha = 0.8) + 
  geom_line(data = p, aes(year, exp(est)), color = pal[2], linewidth = 0.8) +
  geom_ribbon(data = p, aes(x = year,
                            ymin = exp(est - 1.96*est_se),
                            ymax = exp(est + 1.96*est_se)),
              fill = pal[1], alpha = 0.3, inherit.aes = FALSE) +
  facet_grid2(metric~species,
              scales = "free_y", independent = "y") +
  theme(aspect.ratio = 6/7,
        strip.text.x.top = element_text(size = 10, margin = unit(rep(0.06, 4), "cm")),
        strip.text.y.right = element_text(size = 10, margin = unit(rep(0.06, 4), "cm"))) +
  labs(x = "Year", y = "Predation")

tag_facet(p0, fontface = 1, col = "gray30", size = 3)

ggsave(paste0(home, "/figures/pred_yearly.pdf"), width = 17, height = 9, units = "cm")  

Plot weighted predation in space

spatial_pred <- bind_rows(
  get_spatial_predation(sim_her, years = c(1994, 2022), threshold = 1, power = 2, prey_name = "Herring"),
  get_spatial_predation(sim_sad, years = c(1994, 2022), threshold = 1, power = 2, prey_name = "Saduria"),
  get_spatial_predation(sim_spr, years = c(1994, 2022), threshold = 1, power = 3, prey_name = "Sprat")
)

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: no rows removed

mutate: new variable 'pred' (double) with 16,286,659 unique values and 0% NA

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,574 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: removed 3,897 rows (<1%), 16,282,762 rows remaining

mutate: new variable 'pred' (double) with 16,282,762 unique values and 0% NA

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,574 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: removed 26,009 rows (<1%), 16,260,650 rows remaining

mutate: new variable 'pred' (double) with 16,260,650 unique values and 0% NA

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,574 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

plot_map2 <- plot_map +
  guides(fill = guide_colorbar(position = "inside")) +
  theme(legend.position.inside = c(0.23, 0.84),
        legend.key.width = unit(0.2, "cm"),
        legend.key.height = unit(0.25, "cm"),
        legend.text = element_text(size = 6),
        legend.direction = "vertical",
        legend.title = element_text(size = 8),
        plot.title = element_text(margin = margin(0,0,-10,0))) +
  facet_wrap(~year, ncol = 1)

annotations <- data.frame(
  xpos = c(-Inf,-Inf),
  ypos =  c(-Inf, Inf),
  year = c(1994, 2022),
  hjust = c(-0.3, -0.3),
  vjust = c(-20, 1.5))

# Max predation by species (because I trim the plot)
spatial_pred |> 
  summarise(max = max(pred_mean), .by = species)

summarise: now 3 rows and 2 columns, ungrouped

# A tibble: 3 × 2
  species     max
  <chr>     <dbl>
1 Herring    3.19
2 Saduria 1066.  
3 Sprat    239.  

# Plot!

# Herring
p1 <- plot_map2 + 
  geom_raster(data = spatial_pred |> filter(species == "Herring"),
              aes(X*1000, Y*1000, fill = pred_mean)) + 
  geom_sf(color = "gray80") + 
  scale_fill_viridis(trans = "sqrt", name = "Predation\n(kg)",
                     option = "cividis",
                     na.value = "#FFEA46FF",
                     limits = c(0, quantile(filter(spatial_pred, species == "Herring")$pred_mean, 0.999))) +
  labs(title = "Herring") +
  theme(axis.title.x = element_blank()) + 
  geom_text(data = annotations |> mutate(annotateText = c("(a)","(d)")), size = 3.2,
            aes(x = xpos, y = ypos, hjust = hjust, vjust = vjust, label = annotateText))

filter: removed 65,148 rows (67%), 32,574 rows remaining

filter: removed 65,148 rows (67%), 32,574 rows remaining

mutate: new variable 'annotateText' (character) with 2 unique values and 0% NA

# Saduria
p2 <- plot_map2 + 
  geom_raster(data = spatial_pred |> filter(species == "Saduria"),
              aes(X*1000, Y*1000, fill = pred_mean)) + 
  geom_sf(color = "gray80") + 
  scale_fill_viridis(trans = "sqrt", name = "Predation\n(kg)",
                     option = "cividis",
                     na.value = "#FFEA46FF",
                     limits = c(0, quantile(filter(spatial_pred, species == "Saduria")$pred_mean, 0.999))) +
  labs(title = "Saduria") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank()) +
  geom_text(data = annotations |> mutate(annotateText = c("(b)","(e)")), size = 3.2,
            aes(x = xpos, y = ypos, hjust = hjust, vjust = vjust, label = annotateText))

filter: removed 65,148 rows (67%), 32,574 rows remaining

filter: removed 65,148 rows (67%), 32,574 rows remaining

mutate: new variable 'annotateText' (character) with 2 unique values and 0% NA

# Sprat
p3 <- plot_map2 + 
  geom_raster(data = spatial_pred |> filter(species == "Sprat"),
              aes(X*1000, Y*1000, fill = pred_mean)) + 
  geom_sf(color = "gray80") + 
  scale_fill_viridis(trans = "sqrt", name = "Predation\nkg)",
                     option = "cividis",
                     na.value = "#FFEA46FF",
                     limits = c(0, quantile(filter(spatial_pred, species == "Sprat")$pred_mean, 0.999))) +
  labs(title = "Sprat") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        axis.title.x = element_blank()) +
  geom_text(data = annotations |> mutate(annotateText = c("(c)","(f)")), size = 3.2,
            aes(x = xpos, y = ypos, hjust = hjust, vjust = vjust, label = annotateText))

filter: removed 65,148 rows (67%), 32,574 rows remaining

filter: removed 65,148 rows (67%), 32,574 rows remaining

mutate: new variable 'annotateText' (character) with 2 unique values and 0% NA

p1 + p2 + p3 + plot_layout(axes = "collect")

ggsave(paste0(home, "/figures/spatial_predation.pdf"), width = 19, height = 13, units = "cm")

Plot CV in space for cod density and predation

# Apply functions and combine
cv_pred <- bind_rows(
  get_cv_predation(sim_her, years = c(1994, 2022), threshold = 1, power = 2, prey_name = "Herring"),
  get_cv_predation(sim_sad, years = c(1994, 2022), threshold = 1, power = 2, prey_name = "Saduria"),
  get_cv_predation(sim_spr, years = c(1994, 2022), threshold = 1, power = 3, prey_name = "Sprat")
)

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: no rows removed

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,570 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: removed 3,897 rows (<1%), 16,282,762 rows remaining

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,570 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (sim, value) [was 504897x532, now 252448500x34]

filter: removed 236,161,500 rows (94%), 16,287,000 rows remaining

mutate: changed 16,287,000 values (100%) of 'value' (0 new NA)

left_join: added one column (sim_density)

           > rows only in x           341

           > rows only in y  (236,136,596)

           > matched rows      16,286,659

           >                 =============

           > rows total        16,287,000

drop_na: removed 341 rows (<1%), 16,286,659 rows remaining

filter: removed 26,009 rows (<1%), 16,260,650 rows remaining

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,570 unique values and 0% NA

        new variable 'species' (character) with one unique value and 0% NA

# Now do cod density
density_cv <- pred_grid_density |> 
  filter(year %in% c(1994, 2022)) |> 
  summarise(density_mean = mean(sim_density),
            density_sd = sd(sim_density),
            .by = c(year, X, Y)) |>
  mutate(cv = density_sd / density_mean) |> 
  dplyr::select(-density_mean, -density_sd) |> 
  mutate(species = "Cod density")

filter: removed 236,136,596 rows (94%), 16,286,659 rows remaining

summarise: now 32,574 rows and 5 columns, ungrouped

mutate: new variable 'cv' (double) with 32,574 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

# Combine
cv <- bind_rows(cv_pred |> 
                  mutate(species = fct_recode(species,
                                              "Relative herring weight" = "Herring",
                                              "Relative Saduria weight" = "Saduria",
                                              "Relative sprat weight" = "Sprat")),
                density_cv)

mutate: converted 'species' from character to factor (0 new NA)

p1 <- plot_map + 
  geom_raster(data = cv,
              aes(X*1000, Y*1000, fill = cv)) + 
  geom_sf(color = "gray80") + 
  scale_fill_viridis() + 
  facet_grid(year ~ species) +
  guides(fill = guide_colorbar(position = "inside")) + 
  theme(legend.position.inside = c(0.035, 0.82),
        legend.key.width = unit(0.2, "cm"),
        legend.key.height = unit(0.25, "cm"),
        legend.text = element_text(size = 6),
        legend.direction = "vertical",
        legend.title = element_text(size = 8),
        strip.text.x.top = element_text(size = 10, margin = unit(rep(0.06, 4), "cm")),
        plot.title = element_text(margin = margin(0,0,-10,0)))

tag_facet(p1, fontface = 1, col = "gray30")

ggsave(paste0(home, "/figures/supp/cv_density_relative_prey.pdf"), width = 19, height = 10, units = "cm")

Spatial overlap over time

# Calculate overlap indices in space
pred_grid_sum <- pred_grid_density |>
  mutate(cod_spr_ovr = loc_collocspfn(pred = sim_density, prey = biomass_spr),
         cod_her_ovr = loc_collocspfn(pred = sim_density, prey = biomass_her),
         cod_sad_ovr = loc_collocspfn(pred = sim_density, prey = saduria),
         .by = c(year, sim)) |> # Stop here for plotting overlap in space
  summarise(cod_spr_ovr_tot = sum(cod_spr_ovr),
            cod_her_ovr_tot = sum(cod_her_ovr),
            cod_sad_ovr_tot = sum(cod_sad_ovr),
            .by = c(year, sim)) |> 
  # Now calculate mean and sd of annual overlap indices across sims
  summarise(cod_spr_ovr_tot_mean = mean(cod_spr_ovr_tot),
            cod_spr_ovr_tot_sd = sd(cod_spr_ovr_tot),
            cod_her_ovr_tot_mean = mean(cod_her_ovr_tot),
            cod_her_ovr_tot_sd = sd(cod_her_ovr_tot),
            cod_sad_ovr_tot_mean = mean(cod_sad_ovr_tot),
            cod_sad_ovr_tot_sd = sd(cod_sad_ovr_tot),
            .by = year)

clean_ovr <- pred_grid_sum |>
  pivot_longer(c(cod_spr_ovr_tot_mean, cod_her_ovr_tot_mean, cod_sad_ovr_tot_mean,
                 cod_spr_ovr_tot_sd, cod_her_ovr_tot_sd, cod_sad_ovr_tot_sd),
               names_to = "group", values_to = "value") |> 
  mutate(var = ifelse(grepl("_sd", group), "sd", "mean"),
         species = ifelse(grepl("her", group), "Herring", "Sprat"),
         species = ifelse(grepl("sad", group), "Saduria", species)) |> 
  dplyr::select(-group) |> 
  pivot_wider(values_from = value, names_from = var) |> 
  mutate(species = as.factor(species))

# beta-GAM
beta_gam <- sdmTMB(mean ~ species + s(year, by = species),
                   spatial = "off", 
                   family = Beta(link = "logit"),
                   data = clean_ovr)

nd <- expand_grid(species = as.factor(c("Sprat", "Herring", "Saduria")),
                  year = unique(clean_ovr$year))

p <- predict(beta_gam, newdata = nd, se_fit = TRUE)

clean_ovr <- clean_ovr |> 
  mutate(lwr = (mean - sd),
         upr = (mean + sd))

p1 <- ggplot(clean_ovr, aes(year, mean)) + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0, color = "grey40", alpha = 0.7) + 
  geom_point(color = "grey30", alpha = 0.8) + 
  geom_line(data = p, aes(year, boot::inv.logit(est)), color = pal[2], linewidth = 0.9) +
  geom_ribbon(data = p, aes(x = year,
                            ymin = boot::inv.logit(est - 1.96*est_se),
                            ymax = boot::inv.logit(est + 1.96*est_se)),
              fill = pal[1], alpha = 0.3, inherit.aes = FALSE) +
  # facet_wrap(~factor(species, levels = c("Herring", "Sprat", "Saduria")),
  #            ncol = 3) +
  facet_wrap(~species, ncol = 3) + 
  theme(aspect.ratio = 6/7,
        strip.text.x.top = element_text(size = 10, margin = unit(rep(0.06, 4), "cm"))) +
  labs(x = "Year", y = "Spatial overlap index\n(Local index of collocation)")

tag_facet(p1, fontface = 1, col = "gray30")

ggsave(paste0(home, "/figures/annual_overlap.pdf"), width = 17, height = 7, units = "cm")

Correlation between spatial overlap and predation

clean_corr <- clean_pred |> 
  left_join(clean_ovr |> 
              rename(mean_ovr = mean,
                     lwr_ovr = lwr,
                     upr_ovr = upr),
            by = c("year", "species"))

p2 <- ggplot(clean_corr, aes(mean_ovr, mean)) +
  geom_point(color = "gray30", alpha = 0.7) +
  geom_errorbar(aes(ymax = upr, ymin = lwr), color = "gray40",
                alpha = 0.3) +
  geom_errorbar(aes(xmax = upr_ovr, xmin = lwr_ovr), color = "gray40",
                alpha = 0.5) +
  facet_grid2(metric ~ species,
              scales = "free", independent = "all", switch = "y") + 
  theme(aspect.ratio = 1,
        strip.placement = "outside",
        strip.text.x.top = element_text(size = 10, margin = unit(rep(0.04, 4), "cm")),
        strip.text.y.left = element_text(size = 10, margin = unit(rep(0.07, 4), "cm"))) + 
  labs(y = "", x = "Spatial overlap") +
  stat_correlation(size = 2, label.x = 0.98,
                   use_label(c("R", "R.CI")))

tag_facet(p2, fontface = 1, col = "gray30", size = 3)

ggsave(paste0(home, "/figures/corr.pdf"), width = 17, height = 12, units = "cm")

Spatial overlap in space

pred_grid_sp <- pred_grid_density |>
  drop_na(saduria) |> 
  mutate(cod_spr_ovr = loc_collocspfn(pred = sim_density, prey = biomass_spr),
         cod_her_ovr = loc_collocspfn(pred = sim_density, prey = biomass_her),
         cod_sad_ovr = loc_collocspfn(pred = sim_density, prey = saduria),
         .by = c(year, sim)) |> 
  # Sum across sims, by year and spatial location
  summarise(cod_spr_ovr = mean(cod_spr_ovr),
            cod_her_ovr = mean(cod_her_ovr),
            cod_sad_ovr = mean(cod_sad_ovr),
            .by = c(year, X, Y))

# Sprat
plot_map_fc + 
  geom_raster(data = pred_grid_sp,
              aes(X*1000, Y*1000, fill = cod_spr_ovr)) +
  scale_fill_viridis(trans = "sqrt", name = "Cod-Sprat overlap",
                     option = "mako",
                     na.value = "#DEF5E5FF",
                     limits = c(0, quantile(pred_grid_sp$cod_spr_ovr, 0.999))) +
  labs(caption = paste("maximum estimated overlap =", round(max(pred_grid_sp$cod_spr_ovr), digits = 4)))

ggsave(paste0(home, "/figures/supp/spr_overlap_space.pdf"), width = 19, height = 19, units = "cm")

# Herring
plot_map_fc + 
  geom_raster(data = pred_grid_sp,
              aes(X*1000, Y*1000, fill = cod_her_ovr)) +
  scale_fill_viridis(trans = "sqrt", name = "Cod-Herring overlap",
                     option = "mako",
                     na.value = "#DEF5E5FF",
                     limits = c(0, quantile(pred_grid_sp$cod_her_ovr, 0.999))) +
  labs(caption = paste("maximum estimated overlap =", round(max(pred_grid_sp$cod_her_ovr), digits = 4)))

ggsave(paste0(home, "/figures/supp/her_overlap_space.pdf"), width = 19, height = 19, units = "cm")

# Saduria
plot_map_fc + 
  geom_raster(data = pred_grid_sp,
              aes(X*1000, Y*1000, fill = cod_sad_ovr)) +
  scale_fill_viridis(trans = "sqrt", name = "Cod-Saduria overlap",
                     option = "mako",
                     na.value = "#DEF5E5FF", limits = c(0, quantile(pred_grid_sp$cod_sad_ovr, 0.999))) +
  labs(caption = paste("maximum estimated overlap =", round(max(pred_grid_sp$cod_sad_ovr), digits = 4)))

ggsave(paste0(home, "/figures/supp/sad_overlap_space.pdf"), width = 19, height = 19, units = "cm")

# All years pooled
pred_grid_sp_sum <- pred_grid_sp |> 
  filter(year %in% c(1994, 2022)) |> 
  pivot_longer(c(cod_spr_ovr, cod_her_ovr, cod_sad_ovr),
               names_to = "species") |> 
  summarise(mean_overlap = mean(value), 
            .by = c(year, X, Y, species)) |> 
  mutate(species = fct_recode(species,
                              "Sprat" = "cod_spr_ovr",
                              "Herring" = "cod_her_ovr",
                              "Saduria" = "cod_sad_ovr"))
  
p3 <- plot_map + 
  geom_raster(data = pred_grid_sp_sum,
              aes(X*1000, Y*1000, fill = mean_overlap)) +
  geom_sf(color = "gray80") + 
  #facet_wrap(~factor(species, levels = c("Herring", "Sprat", "Saduria"))) +
  facet_grid(year~species) +
  scale_fill_viridis(trans = "sqrt", name = "Spatial overlap", option = "mako") +
  theme(legend.position.inside = c(0.08, 0.845),
        legend.key.width = unit(0.2, "cm"),
        legend.key.height = unit(0.33, "cm"),
        legend.text = element_text(size = 7),
        legend.direction = "vertical",
        strip.text.x.top = element_text(size = 10, margin = unit(rep(0.06, 4), "cm"))) + 
  guides(fill = guide_colorbar(position = "inside")) +
  NULL

tag_facet(p3, fontface = 1, col = "gray30")

ggsave(paste0(home, "/figures/overlap_space.pdf"), width = 19, height = 14, units = "cm")